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ABSTRACT 


The  manual  method  of  minimum-time  ship  routing  has 
proved  so  successful  over  the  last  several  years  that 
ever  increasing  numbers  of  ships'  captains  are  requesting 
this  service.  The  desirability  of  the  electronic  computer 
for  the  computation  is  therefore  readily  apparent.  In 
this  investigation  the  minimum-time  route  is  determined  by 
calculus  of  variations ,  rather  than  the  conventional  manual 
technique.  The  present  method  consists  of  solving  the 
associated  Euler  equation  by  numerical  integration  on  the 
Control  Data  Corporation  model  1604  digital  computer.  In 
case  1  the  ship's  speed  is  assumed  to  be  primarily  a  func¬ 
tion  of  position;  and  in  case  2  the  ship's  speed  is  taken  to 
be  a  function  of  its  direction  as  well  as  position.  Case 
2  is  recommended  for  operational  adaptation. 

The  writer  wishes  to  express  his  appreciation  for  the 
assistance  given  him  by  Professor  George  J,  Haltiner  of  the 
United  otates  Naval  Postgraduate  ochool  in  this  study. 
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T-J3LE  07  JiriBCLo 


T  =  Time,  hours 

V  =  Chip's  speed,  knots 

4'  =  l/V,  hours  per  nautical  mile 

ds  =  Element  of  arc  length  along  the  track 
y '  =  dy/dx 

y"  =  d2y/dx2 

y^  =  y  value  at  the  itt  grid  point 

y.+j  =  y  value  at  the  (i+J)  grid  point 

7j__l  =  y  value  at  the  (i-l)  grid  point 

D  =  Finite  increment  in  the  x-direction 

=  "  i  +  l  “  J  i- J 

2D 


y 


t  t 


n+i  -2n  +yi- 


D 


=  Uave  height,  in  feet 


P  =  Direction  from  which  waves  are  coming 


<  =  Heading  of  ship 


=  Fetation  angle  between  cartesian  coordinates 
of  wave  model  (x,r,  y  )  and  of  ship 
routing  (x,y) 
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Introduction 


j. . 

.'i.  manual  method  of  minimum-time  ship -routing  was 
developed  and  tested  by  the  U.b.  Navy  Hydrographic  Office 
[l] .  In  the  investigation  it  was  discovered  that  the  most 
important  parameter  in  retarding  a  ship's  progress  was 
wave  action.  This  action  manifests  itself  in  two  ways; 
first,  the  direct  reduction  by  the  resistance  of  the  water 
to  the  ship's  passage,  and  second,  the  indirect  reduction 
through  the  voluntary  decrease  of  ship's  speed  to  reduce 
the  violent  effects  of  roll,  pitch  and  heave.  The  two 
most  important  components  of  this  wave  action  are  wave 
height  and  wave  direction  t2]  . 

The  Fleet  Weather  Facilities  at  .'.lameda  and  at  Norfolk 
have  been  very  successful  with  the  operational  application 
of  the  manual  method  to  hundreds  of  cases  during  the  last 
several  years  [l]. 

With  the  advent  and  continued  improvement  of  elect¬ 
ronic  computers,  there  has  been  a  simultaneous  increase  in 
their  application  to  scientific  and  engineering  problems. 
Cne  such  example  is  the  processing  of  meteorological  data 
and  the  subsequent  forecasts. 

The  application  of  the  electronic  computer  to  the 
minimum-time  ship  routing  problem  appears  desirable  as 
ever-increasing  numbers  of  ships'  captains  are  requesting 
this  service.  The  first  natural  _ impulse  would  be  to  adapt 
the  present  manual  method  to  the  computer.  However  there 
is  a  branch  of  mathematics  called  calculus  of  variations 
which  deals  with  the  problem  of  minimizing  or  maximizing 
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the  value  of  an  integral  by  the  determination  of  an  approp¬ 
riate  function,  in  this  case  the  proper  route  which  will 
minimize  the  time  required  for  a  ship  to  travel  between  two 
ports  of  call. 

Under  quite  general  conditiors  it  may  be  shown  that 
the  function  y(x)  necessary  to  minimize  or  maximize  the 


integral 


r*2 

J  F (yyyjy/)  dx 


CD 


must  fulfill  the  Euler  equation 

d-  M  -  M  =  o 

dx  £y' 

The  boundary  conditions  applied  here  are  that  the  y  values 
at  x,  and  xe  are  fixed  [3]. 

The  time  T  required  for  a  ship  or  other  vehicle  to 
transverse  a  track  is  given  by  the  line  integral 


(2) 


T  = 


ds 


In  cartesian  coordinates  ds  may  be  replaced  by  [dx"  + 
and  (1)  may  now  be  put  in  the  form 

+(y'f]/k  Jx 


(3) 


'.mere :  f  [  I  +  (y')*]  =  p (X)  y}  y')  . 

The  Euler  equation  (2),  which  must  be  satisfied  for  a 
minimum,  becomes 

j  ilHi+y*)*]  _  Mni  + ?’*)*}  =  0 

dx  <$y'  &  y 


(4) 


(5) 
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.Several  methods  are  available  for  the  solution  of  the 
huler  equation.  One  possibility  is  to  treat  it  as  on 
initial- value  problem.  For  a  given  initial  heading  of  the 
ship  there  w ill  be  an  error  in  the  arrival  point  as  det¬ 
ermined  by  this  method  from  the  desired  destination.  Thus 
a  first  estimate  is  made  for  the  initial  heading  and  the 
subsequent  error  is  determined.  A  second  initial  heading 
is  then  chosen  in  such  a  way  as  to  reduce  the  corresponding 
error  at  the  final  point.  Thus  the  initial  heading  is 
successively  modified  until  the  destination  error  is  within 
acceptable  limits  of  accuracy. 

A  relaxation  method  is  used  in  this  study  and  will  be 
described  in  a  later  section. 

Si chard  W .  J ame  s  [2]  found  that  theoretical  equations 
representing  ships'  speed  under  various  seas  did  not  fit 
observed  data  taken  from  ships’  logs.  Cne  of  the  many 
reasons  given  for  this  discrepancy  was  the  human  factor  of 
voluntarily  reducing  the  ship's  speed  in  order  to  reduce 
the  motions  of  the  ship.  James  first  classified  his  data 
according  to  four  90-degree  sectors  of  relative  wave  dir¬ 
ection.  Waves  that  approached  a  ship  within  45  degrees  of: 

1.  the  bow  were  classified  as  head  waves; 

2.  the  beam  were  classified  as  beam  waves; 

3.  the  stern  were  classified  as  following  waves. 

Then  winder  each  classification,  wave  height  (feet)  versus 
ship's  speed  (knots)  was  plotted.  Using  the  method  of  least 
squares  on  the  assumed  straight  line  function  V=  a  +  bH, 
the  constants  "a"  and  "b"  were  determined  to  give  tlie  best 


fit  for  each  of  the  three  classifications. 


The  resulting 


equations  are: 

Head  Geas  V  =  18.29  -  0.529  H 
Bcara  8c as  V  =  17.90  -  0.276  H 
Follov/ing  3eas  7  =  18.50  -  0.168  H 


(6) 


Note  that  zero  wave  height  does  not  result  in  the  same 


speed  in  the  three  equations. 
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3. 


■) 

I  ,  case  1  it  is  assunec  t.  at  t’  e  shi_  9  s  s^eccl  jl  a 
function  of  _.osition  only  and  not  a  function  of  ship's  head¬ 
ing.  Strictly,  this  implies  ship's  speed  is  a  function  of 
wave  height  alone,  however,  in  practice  the  empirical 
formulas  (6)  give  no  variation  in  speed  for  a  considerable 
change  of  ship's  heading  in  their  respective  broad  sectors 
of  head,  beam,  and  following  seas.  Hence  the  ship's  dir¬ 
ection  may  be  changed  appreciably  without  altering  the 
empirically-determined  ship's  speed.  It  follows  that  in 
many  cases  no  great  error  would  arise  from  neglecting  those 
ter  is  in  the  iuler  equation  involving  the  variation  of  the 
shirk's  speed  as  the  course  changes  slightly  in  the  success¬ 
ive  approximations,  nevertheless  this  more  general  case  is 
considered  in  the  next  section. 

If  V,  and  therefore  i<  ,  is  a  function  of  x  and  y  but  not 
y' ,  the  Juler  equation  (5)  becomes 

f  y"  +  (ix  r-  +  y,z)  =  o  (?) 

^  A 

Using  central  finite  difference  approximations  y'  and  y' ' 
for  y'  and  y* '  in  (7);  the  equation  for  y^  becomes 


dince  the  ship's  speed  is  expressed  analytically  (6), 
the  derivatives  may  be  expressed  analytically  and  evaluated; 
or  the  derivatives  may  be  approximated  by  finite  differences. 
The  latter  method  wo  Id  probably  be  usee1  in  practical 
operations. 
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Equation  (8)  is  solved  by  the  relaxation  method,  which 
consists  of  calculating  the  nb  approximation  for  y  ,  denoted 

(  n  ^ 

as  y^v  '  ,  using  the  previously  calculated  approximations 

y.  and  y.  The  initial  estimate  is  somewhat 

arbitrary  and  nay  be  taken  as  the  great  circle  route,  for 
example.  The  initial  and  final  points  remain  as  fixed  bound¬ 
ary  conditions. 

International  steamer  tracks,  which  are  modifications 
of  routes  laid  out  by  Lieutenant  Ilatthew  F.  Taury,  U32* 

(circa  1850),  are  quite  often  selected  by  ships'  captains 
when  the  minimal-time  ship  routing  service  is  not  used. 

These  routes  were  the  best  statistical  routes,  based  on  the 
climatology  of  the  storm  tracks  and  sea  conditions,  to 
revent  damage  and  to  decrease  sailing  time  [1] . 

However,  to  keep  the  problem v  general,  the  initial 
estimate  for  the  correct  route  is  taken  as  the  great  circle 
which  passes  through  the  point  of  departure  and  the  point 
of  arrival.  To  take  full  advantage  of  this  choice,  a 
Lambert  Conformal  chart  is  used  on  which  a  great  circle  is 
approximately  a  straight  line  [4] .  The  secant  projection 
of  the  Lambert  Conformal  chart  is  chosen  in  order  to  have  a 
broader  area  of  minimum  distortion.  The  Image  scale  e  , 

which  is  the  dimensionless  ratio  of  image  distance  to  earth 

/ 

distance  is  given  by  .  ‘ 


£ 


Sin  <P, 

Sin4> 


-ta.n  -jr 
tan  & 


n 
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n  Jjj-g-Sm  4\  JLrrp  Sm  $2 

./here:  =  — j-* - -r - a - — 

With  the  standard  parallels  of  j,  =  30  degrees  and  ^  =  60 

degrees,  n  =  0.716  [53*  Then  e  =  1.0  at  30  and  60  degrees. 
The  minimum  value  of  e  is  at  45  degrees  and  is  approximately 
0.97.  By  staying  within  a  range  of  about  25  to  65  degrees 
of  latitude,  e  may  be  neglected  with  an  error  of  3 %  or  less. 
The  map  scale  does  not  enter  into  the  computation  of  the 
problem. 

In  order  to  test  the  method,  a  fairly  general  hypo¬ 
thetical  wave  distribution  was  designed,  as  shown  schematic¬ 
ally  in  figure  1.  In  the  central  region  between  the 
dividing  lines  W1  and  W2,  the  wave  height  is  a  function  of 
yw  together  with  the  maximum  wave  height  which  is  constant 
along  the  dividing  line  ¥3.  Cutside  this  central  area  the 
wave  height  is  a  function  of  x  and  y  ,  with  wave  height  de- 
creasing  aw a y  from  the  central  region.  To  avoid  discon¬ 
tinuities  narrow  linear  transition  zones  of  speed  were 
inserted.  The  computational  procedure  is  programmed  in  such 
a  way  that  the  various  parameters  defining  the  model  may  be 
changed  at  will.  Thus  a  large  variety  of  sea  conditions 
may  be  simulated. 
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_  _  w3  BEAM  SEA 


CASE  1  WAVE  MODEL 


Figure  1 


3*  Case  2. 

I  .  t  is  section  the  more  general  case  where  the  shi  's 
si ecd  varies  continuously  with  relative  wave  direction  will 
be  considered.  Thus  i  ~  ,{'(x,y,y')  ,  and  the  Euler  equation 
(5)  reduces  to 


Ml 

bx 


y' 


+ 


< at 

by 


ML 

byj 


y'  + 


by'bx 

by* 


o 


(9) 


When  this  equation  (9)  is  expressed  in  finite  differ¬ 


ence  form  and  solved  for  y^  it  becomes 


+ 


ax  y  '  ~5y  | 

;  . 

[ft  , 

[bybx  by 'by  ^ 

oo) 

U[j+P/2]“«+2 

$,+  r 

by' '  by'zl 

i+y*] 

h 

In  order  to  allow  for  a  continuous  change  of  ship's 
speed  with  relative  wave  direction,  James'  three  equations 
(6)  may  be  approximated  by 


V  =  18.23  -  0.302  +  0.189  Cos  (p-cQ  H  C11) 

If  a  better  approximation  with  head  and  beam  wave  conditions 
is  desired,  with  a  slight  loss  of  accuracy  with  following 
wave  conditions,  then  (11)  may  be  modified  to 

V  =  18.23  -  0.302  +  0.225  Cos  (p-oQ  H  (12a) 

On  the  other  hand,  a  better  approximation  with  following 

and  bean  wave  conditions,  with  a  slight  loss  of  accuracy 
with  head  \;ave  conditions,  is 
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V  =  18.23  -  0.302  +  0.153  Cos  (f.-oO  H  (12b) 

TI-J  last  equation  was  used  for  the  results  in  table  1; 
however  the  three  constants  in  order  were  rounded  to  18.0, 
0.30  and  0.15,  respectively,  bTith  these  values,  the  three 
cases  of  head,  bean  and  following  seas  are  represented  by 

Head  Teas  7=  18.0  -  0.45  H 

Beara  Seas  V=  18.0  -  0.30  II  (13) 

Following Seas  7=  18.0  -  0.15  H 

In  order  to  tost  the  effect  of  including  the  variation 
of  ship's  speed  with  its  direction  in  the  duler  equation, 
only  the  addition  of  wave  direction  at  every  point  to  the 
previously-designed  wave  model  was  necessary,  as  shown 
schematically  in  figure  2.  Outside  the  central  region  of 
the  odel  the  wave  direction  follows  a  circular  cl  uige  until 
it  returns  to  the  central  area.  ..  transition  zone  is  needed 
only  to  separate  head  and  following  seas.  This  is  accomp¬ 
lished  by  a  rapid  change  of  wave  direction  between  these 
two  areas.  The  remaining  short  line  of  discontinuity  of 
speed  is  positioned  at  the  start  of  each  specific  problem 
near  the  center  of  the  transition  zone,  but  in  such  a 
position  that  it  will  not  become  involved  in  the  subsequent 
calculations. 

Other  conditions  remain  as  before. 
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CASE  2  WAVS  MODEL 


i,  xeswlts. 

Case  1  :n -1  c:  2  results  were  ’oof.  obtui.  .ed  udth 

initial  conditions  as  s  ccified  beloi  (suo  figures  1  and  2), 
The  great-circle  distance  between  dep.  rt  are  end  arrival 
points  was  4,800  n  mi.  The  finite  increment  in  tiie  re¬ 
direction  D  was  60  n  ni,  She  rotation  angle  &-t was  zero; 
thus  the  two  coordinate  axes  coincided.  was  1200  and  Wg 

was  3600  n  mi.  The  general  wave-directional  pattern  was 
cyclonic,  that  is,  head  seas  above  the  great  circle  (pos¬ 
itive  y-direction)  and  following  seas  below  the  great 
circle  (negative  y-direction).  nil  initial  y  values  were 
taien  as  zero  along  the  great  circle;  therefore  the  maximum 
displacements  shown  in  the  table  are  the  perpendicular 
distances  away  ^ ron  the  great-circle  route.  The  maximum 
wave  height  was  varied  from  zero  to  thirty  feet  at  one-foot 
increments,  have  height  decreased  at  the  rate  of  0.03  feet 
per  n  mi  with  respect  to  x  and  0.06  feet  per  n  mi  with 
respect  to  y. 

The  calculated  results  of  cases  1  and  2  varied  smoothly, 
therefore  only  the  results  obtained  at  five-foot  increments 
of  maximum  wave  height  are  listed  in  table  1.  The  results 
in  the  table  are  rounded  to  the  nearest  one-tenth  unit. 

In  case  1  the  great-circle  estimate  went  through  the 
middle  of  the  head-following  linear  transition  zone,  that 
is,  Wj  was  zero  n  mi.  The  apparent  discrepancy  for  zero 
w_.ve-height  is  a  result  of  the  unequal  speeds  from  the 
James'  equation  (6)  for  zero  wave  height,  and  the  passage 
through  a  linear  transition  zene  of  the  model. 
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In  case  t'  o  re  t-circle  estimate  went  through,  on 
t.  v  u  _  or  li  .  it  01  t  o  tr •  j.  sition  zone,  that  is,  head  seas 
were  encountered  throughout  the  central  region,  with  set 
at  ninus  one  n  mi. 

In  comparing  the  results  of  cases  1  and  2,  the  minimum 
times  and.  maximum  displacements  should  he  compared  rather 
than  the  great- circle  times,  since  the  greatest  differences 
between  equations  (6)  and  (12b)  or  (13)  occur  with  head  seas. 

In  cases  1  and  2  the  convergence  to  the  minimum  time 
varied  in  a  nonlinear  manner  with  the  rate  of  convergence 
decreasing  as  the  minimum-time  track  was  approached.. 

To  make  certain  that  the  calculated  results  were  not 
just  relative  minimum s  due  to  some  unforeseen  combination 
of  speed  equations  and  wave  model,  initial  estimates  other 
than  the  great  circle  were  used.  These  vjere  taken  on  both 
sides  of  the  computed  mininum-time  track;  nevertheless  the 
final  track  was  the  sane. 
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RESULTS  OP  CASE  1  AND  CASE  2 


C-SS73 

1 

WAVE 

CASE  2 

TILE  (HOURS) 

TIALIIIUM 

LEIGLT 

l alii  run 

TILE  (HOURS) 

GRS-T 

CIRCLE 

KILIiRJU 

ROUTE 

DIOPLaCSLENT 
(  N.  Liles)) 

(FEET) 

DISPLACEL^.'T 
(  IT.  Liles  ) 

IIIITILUL 

ROUTE 

GREAT 

CIRCLE 

264.8 

265.6 

-0.7 

0 

0.0 

286.7 

266.7 

287-8 

260 . 2 

-114.8 

5 

-115.3 

267.1 

286.4 

321.7 

265-0 

-196.5 

10 

-196.2 

288.0 

314.7 

373-5 

266.5 

-280.7 

15 

-276.8 

269-5 

355.4 

460 . 3 

268.6 

-361.1 

20 

-361.8 

271.5 

416.9 

632.8 

271.2 

-443 . 2 

25 

-443 . 2 

274.0 

517.8 

1156.7 

270.5 

-544.7 

30 

-527.5 

277.2 

713.0 

Table  1 


-v  j-.*.  C  _i~  o  -L  u.  ■  e 

w  calcvlu  j  ' 

1 1  - 


ianxi  ^  ^  jf.v  oo  .minuii- 
n.ius  Co  u  .  "otl .1  i:  min 


was  'i  "covered  th:  t 
xae t hoc' . 

In  many  cm_i..  ccrf.  3  and  set. 


IS  SOUny 

ould  bar  o  eraticnal  ad  a  tatio. :.  of  this 


.c  rollers  solved  by 


relaxation,  it  is  agreed  that  the  best  approximation  to  the 
correct  answer  is  the  re  an  betvreen  the  values  obtained  by 
relaxing  fro:1,  initial  estimates  on  bot.  sides  of  the  cal¬ 
culated  answers.  'owever  in  this  investigation  the  first 
calculated  answer  from  the  great  circle  route  is  well 
within  present-day  ne.vi rational  standards,  and.  no  furt'  cr 
calculation  appears  necessary. 

In  this  study,  case  2  took  about  four  tines  as  long  to 
calculate  as  case  1,  but  since  the  solution  time  involved 
was  still  only  about  to:  minutes  on  the  Control  Data  Corp¬ 
oration  model  1504-  digital  computer,  the  added  generality 
seeraed  worth  the  additional  time.  .Iso  since  ship’s  head¬ 
ing  and  wave  direction  must  be  calculated  in  the  o~  erational 
use  of  case  1  in  order  to  determine  the  roper  speed ,  Vis 
information  may  as  well  be  used  to  its  fullest  extent  by 
utilizing  case  2. 

Che  tire  variation  of  the  sea  may  be  included  in  case 
2  in  much  the  sane  manner  as  the  directional  variation  of 
the  sea  was  to  c^se  1.  fortunately  the  minimum-time  route 
is  generally  away  from  active  storm  centers  and  along  the 
fringe  areas  where  tv  e  actual  thic  variation  of  the  seas  is 


jreatly  reduced .  ince  sea  conditions  decay  slowly,  but 
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nay  build  fairly  ra  icily, 


tixc  er  ror  co  Id  "be  of  or  c 


si^mific;  ce  nor  r  ra_ idly-developing  ctorn  centers. 


Thus  the  only  real  limitation  to  the  a;  plic 
this  method  ay years  to  he  in  the  cuality  of  the 
sec  conditions. 


tion  of 


orecast  of 


Janes'  speed  equations,  anc  their  approximations ,  uere 
used  in  this  study  only  to  provide  realistic  values,  and 
in  an  operational  application  each  ship  could  have  its  own 
characteristic  speed  epuation  inserted  into  the  program. 

JLe  ~ amber t  Conformal  chart,  while  convenient  for  this 
experiment,  is  hy  no  means  necessary,  end  an  operational 
program  should  be  chosen  to  use  the  data  available  from 
forecast  wave  height  and  clirection  most  economically. 


,.s  mentioned  earlier,  a  further  refinement  of  this 
method,  of  minimum-time  ship  routing  would  be  the  inclusion 
of  the  terns  in  the  fuler  equation  that  involve  the  chan'e 
of  wave  height  and  direction  with  time.  jince  the  results 
were  changed  very  little  by  the  additi  ,m  of  the  relative 
wave-direction  terms  in  the  duler  epuation,  in  most  sit¬ 
uations  there  probably  would  be  little  change  with  the 
addition  of  the  time -varied; ion  terms,  even  though  they  are 
quite  numerous.  The  ad.ded  time  required  for  solution,  as 
well  as  the  extra  storage  of  sea-condition  data  in  or  out¬ 
side  of  the  computer  and  the  call  end  recall  of  the  data 
for  at  least  three  time  periods  associated,  with  each 
calculation  of  y^ ,  nay  prove  to  be  prohibitive  from  an 
a  ©ration ° 1  stand  oint. 
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